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This is the Final Technical Report on research conducted under NASA Grant NAGW-1077 
entitled "Mars Tectonics and Volcanism." This grant supported the research of one Investigator 
(S.C. Solomon) and one Ph.D. student (P.J. McGovern) on behalf of the NASA MEVTV (Mars: 
Evolution of Volcanism, Tectonics, and Volatiles) Project over the three-year period 1 June 1987- 
31 May 1990. 

The focus of the research was on three broad areas: 

1) The relation between lithospheric stress in the vicinity of a growing volcano and the 
evolution of eruption characteristics and tectonic faulting; 

2) The relation between elastic lithosphere thickness and thermal structure; 

3) A synthesis of constraints on heat flow and internal dynamics on Mars. 

The best means of conveying the principal findings of our research is through the publications 
reporting our work. On the following page is a list of publications supported under this project 
during the past 3 years. Included in the list are 3 papers, 2 book chapters, 7 abstracts of oral 
presentations, and 1 technical report. 

The remainder of this report consists of two of the most recently completed papers supported 
by this project. The other publications listed are or have been included in other progress reports. 
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Derived values of the thickness the 

estimates of lithospheric thermal gradients and surf ^ ubject to assu med strain rates 

elastic-plastic plate having the same bend g mantle mate rial. Local thermal gradients and 

and temperature-dependent flow laws for cru d m -i f f ctively> at the time of 

heat flow values so estimated were 10-14 K km ana zp > A lba Patera, while 

formation of flexurally induced « r ff n K km ‘ and .7-24 mW m'f, respectively, 
gradients and heat flow values of es mascon an d Olympus Mons at the time of emplace- 

characterized the lithosphere beneath the Isidis m k y f , h elastic lithosphere inferred to 

ment of these loads. On the basis ^^le heafproducfton obtained from SNC meteorites, it is 

support the Tharsis rise and est, mates f^ntle heat productro^ ^ ^ ^ ^ , 5 _ 25 mW m 

suggested that the present average globeil heat ch has been contributed by excess 

Approximately 3-5% of this heat ux volcanic provinces. Most likely, this excess heat flux 

conducted heat in the central regions of mantle plumes. The fractional mantle hea 

has been delivered to the base of he Ubmpbrn by P ^ sjmjlar to that at present 

transport contributed by plumes during the last - o.y. on 
on Earth. 


Introduction 

The thickness of the elastic lithosphere on a planet is 
isentially a measure of the reciprocal of the vertical therma 
■adient in the lithosphere, i.e., the depth to a temperature a 
hich ductile behavior replaces brittle behavior at typi 
eological strain rates. Under flexure there is an elastic 
core” of the lithosphere occupying the depth interval over 
/hich the bending stress is less than an envelope of 
‘strength” versus depth defined by a frictional failure curve 
t shallow depths and a ductile 

Goetze and Evans, 1979; Brace and Kohlstedt, 1980]. At the 
hallowest depths, lithospheric bending leads to faulting to a 
lepth that is dependent on the load, the flexural ngidUy.tmd 
he failure law. The depth of the lower limit to elastic 
yehavior is governed primarily by temperature and also y 
strain rate, Composition, and load magnitude. Estimates of 
elastic lithosphere thickness derived from simple models of 
flexure have been quantitatively related to the verticaly 
averaged thermal gradient of the lithosphere on the Earth 
[e.g., Caldwell and Turcotte, 1979; McNutt, .1984; McA *°° 
et al 1985; Willett et al., 1985; Kusznir and Earner, 1985] 
and Moon [Solomon, 1985], and similar concepts have been 
used to constrain the thickness of the elastic lithosphere o 
Venus [e.g., Solomon and Head, 1984]. In this paper we 

apply these concepts to Mars. .. 

We begin with a review of estimates of the effective 
thickness of the Martian elastic lithosphere. We then convert 
these thickness values to estimates of lithospheric therma 
gradients and heat flow by means of temperature-depend 
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strength envelopes. On the basis of the locations and geo- 
logical epochs appropriate to each estimate of thermal gra- 
£ ^ delate the results to global heat flux, interior 
thermal evolution, lithospheric reheating mechanisms, and 
the evolution of major volcanic provinces on Mars. 

Elastic Lithosphere Thickness 
The thickness T e of the elastic lithosphere of Mars has 

been estimated from the tectonic response to tndiv^dual 

loads [Thurber and Toksdz, 1978; Comer et al., W5 Janie 
and Jannsen, 1986] and from the global response to 
long-wavelength load of the Tharsis rise [Willemann and 

tIZ, mi, « al. m2; Sleep and PM,p , 

1985], A summary of these results is given ir .Table E It is 
important to distinguish the elastic l>thosphere from^e 
thermal or compositional lithosphere inferable from Pratt 
isostatic compensation models [e.g.. Sleep and Phi tops, 
1979] or from the heights of volcanic constructs [V g , > 

Carr 1976]. The base of the elastic or mechamcal l'tho 
sphere is governed by the temperature at which ductile 
strength becomes less than some threshold value at geolog- 
ical strain rates [Goetze and Evans, 1979; Brace and Kohl- 

V,< The 'radial distances of graben circumferential to the 
major volcanoes Ascraeus Mons, Pavon.s Mons, Ars a 
Mons, Alba Patera, and Elysium Mons (Figure l)m >«• *te 
best fitting values for flexural rigidity D of 10 -10 N m at 
the times of graben formation [Comer et al., 1985], these 

are equivalent (for values of the ^hosphenc 

Young’s modulus and Poisson s ratio of E - 1« ^ GPa i an 
„ = 0.25, respectively) to T e in the range i^o krn ln 
addition to best fitting values, Comer et al . [1985] obta 
lower and upper bounds to T„ from formal error analysis 
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iponding to at least the late Amazonian and P erha P s t0 
iomewhat earlier epochs as well [Scott and Tanaka, 1986, 

Tanaka et al, 1988]. Th - 

A variety of models for the deep structure of the Tharsis 
rise have been proposed [e.g.. Sleep and Phillips 979, 
1985; Banerdt et al, 1982; Willemann and Turcotte, 1982 
Finnertv et al, 1988]. The distribution and orientation ot at 
least some of the tectonic features of the Tharsis province 
are best explained if the long-wavelength topography of the 
Tharsis rise was supported for a significant time y t e me 
strength of the global elastic lithosphere 
Turcotte, 1982; Banerdt et al., 1982; Sleep and Phillip , 
19851 Quantitative models of this support that provide a 
reasonable fit to the geoid and to the distribution of tectonic 
features have an elastic lithosphere 100 to 400 km thick 
corresponding to D = 10 25 to 7 x 10 26 N m [Wdlemann and 
Turcotte, 1982; Banerdt et al., mi). The t.me at which this 
estimate is an appropriate average for the planet depends on 
the evolution of the mechanism of support of long 
wavelength topography [Sleep and Phillips, 198 ] an can 
range from Middle Noachian (3.8 b.y. ago) to Upper Am - 
zonian [Tanaka, 1986]. 

The best fitting values for D derived for various loads span 
at least 2 orders of magnitude (Figure 2), and the correspond- 
ing values for T e span nearly a factor of 10 (Figure 3). As is 
apparent from Table 1 and from the discussion above, these 
values are not consistent with a simple progressive increase 
with time in the thickness of the elastic lithosphere of Mars 
The largest estimates of T e , for instance are for the oldest 
(Isidis mascon) and youngest (Olympus Mons) ocal i ho- 
spheric loads considered. Spatial variations in elast.c litho- 
sphere thickness must therefore have been at 'east as impor- 
tant as temporal variations [Comer et al., 198 ]. 
particular, there appears to have been a dichotomy in 
lithosphere thickness that spanned a significant interval of 
time with comparatively thin elastic lithosphere (T e 
km) beneath the central regions of major volcanic provinces 
and substantially thicker elastic lithosphere (T e in excess o 
100 km) beneath regions more distant from volcanic prov- 
ince centers and appropriate for the planet as a whole. 


Thermal Gradients 

Procedure 

All of the values of T e described above were obtained 
under the assumption that the elastic lithosphere of Mars 
behaves as a uniform elastic plate or shell. A etter mo > 
for the lithosphere is an elastic-plastic plate, 
strength as a function of strain rate, platc^curva urc a 
depth [Goetze and Evans, 1979; Brace and ^MsUdtimY 
To estimate the mean lithospheric thermal gradient from 
such an effective elastic lithosphere thickness, T must first 
be converted to the depth T m to the rheological boundary 
marking the base of the mechanical lithosphere (Figure ). 
This conversion is accomplished by adopting a representa- 
tive strain rate and a flow law for ductile deformation of 
material in the lower lithosphere, constructing models of 
bending stress consistent with the adopted streng env 
(ope and finding for each model the equivalent elastic plate 
model having the same bending moment and curvature . . 

procedure has been described by McNutt [1984]. 

The temperature distribution in the mechanical litho- 
sphere is assumed to be locally given by a surface tempera- 
ture 7 of 220 K [Kieffer et al., 1977; Fanale et al., 1982] and 
a uniform vertical gradient dTIdz. We take the representative 
strain rate for the flexural response to each local load I to b 
the quotient of the maximum horizontal strain given by t 

mod®' and t* '<»? °"J dc S 

Of Stratigraphic relations among the pnncipal volcanic 
associated with the major lithospheric loads . [St o t and 
Tanaka, 1986; Tanaka, 1986] and the crater chronology of 
Hartmann et al. [1981], we take the load growth time to be 
10 8 years. The brittle and ductile portions of the strengt 
envelope are approximated by straight lines, so t a e 
bending moment may be found analytically giver ^curva- 
ture and dTIdz [McNutt and Menard, 1982 M. K. McNut , 
personal communication, 1988], The brittle portion of the 
strength envelope is taken from the low-pressure faction law 
of Bverlee [1978], with a lithostatic pressure gradient appro- 
priate to the Martian crust (1 1 MPa km - ’) under the assump- 
tion of negligible pore pressure. If the mechanical litho- 
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effective Young’s modulus is not well known for Mars, to 
first order the resulting relative uncertainty in T e is only one 
third that in E. Tn general, of course, curvature is not a 
constant but varies along a flexural profile; it may be argued 
[ McNutt , 1984], however, that the estimates of D and T e in 
Table 1 are dominated by the portions of the profiles at 
which bending stress, and thus curvature, are near their 
maximum values, so we assume such maximum values in 
using Figures 5 and 6. 

The effective surface temperature T s has an uncertainty of 
perhaps 20-30 K. The mean microwave brightness temper- 
ature of Mars is about 200 K \ Morris on et at 1969]. Climate 
models for Mars predict that the present seasonally averaged 
surface temperature should vary with latitude and should be 
200-230 K at the latitudes (approximately 40°N to 10°S) of 
the features listed in Table 1 [Kieffer et al ., 1977; Female et 
al , 1982]. Such models do not generally include variations in 
topography, which likely would add a superposed variation 
in surface temperature of perhaps 10-20 K [Seiff and Kirk , 
1977] for the range in regional elevations of prominent 
lithospheric loads. An additional uncertainty is contributed 
by poorly constrained differences between the present cli- 
mate of Mars and that at the time of load-induced graben 
formation. A final factor is the possible presence of a 
low-conductivity regolith layer, which can be accommo- 
dated in the simple model adopted here by use of a small 
increase in T s [McNutt and Menard , 1982]. A 20-30 K 
uncertainty in T s contributes only a 3-7% uncertainty to 
estimates of lithospheric thermal gradient and heat flow. 

The assumption of a linear thermal gradient can also be 
questioned. In effect, for a given strain rate and flow law the 
quantity T m essentially gives the depth to a particular 
temperature, and any temperature distribution that passes 
through that temperature-depth point and the surface tem- 
perature T s (and that does not exceed the limiting tempera- 
ture at depths less than T m ) should be regarded as possible. 
Effects such as upward concentration of crustal heat- 
producing elements and a thermal conductivity that in- 
creases with depth will tend to give a temperature distribu- 
tion that is concave downward, so that near-surface 
gradients will be higher than those indicated in Figures 5 and 
6. Further, there is also at least an order or magnitude 
uncertainty in the growth time of the lithospheric loads 
utilized to estimate T e and thus in the effective flexural strain 
rate. A factor of 10 uncertainty in growth time, however, 
contributes only a 4% uncertainty to the derived values of 
T m or dTIdz for either of the adopted flow laws. 

A source of uncertainty is the abundance of water in the 
Martian interior. In the upper crust a significant pore pres- 
sure can lower the frictional strength [Brace and Kohlstedt , 
1980], which would result in modest increases in T m (or 
decreases in dTIdz ) for given values of elastic thickness and 
curvature. Whether water is present in the lower crust and 
upper mantle of Mars is an issue for the estimation of ductile 
strength, since the presence of water can significantly affect 
creep rates in rock [e.g., Ashby and Verrall, 1977]. The 
adopted olivine flow law, although derived from laboratory 
measurements, has been reasonably well validated for the 
terrestrial oceanic mantle and geological strains rates, by 
means of both flexural estimates of mechanical lithosphere 
thickness [ McNutt and Menard , 1982; McNutt , 1984] and 
centroid depths of intra plate earthquakes [Bergman and 
Solomon , 1984], although McNutt and Menard [1982] have 


suggested that the activation energy for creep in the mantle 
may be slightly less than that obtained from laboratory 
measurements on dry olivine. For Mars we note that signif- 
icant water in the mantle and an activation energy less than 
that assumed here [ Goetze , 1978] would both result in lower 
estimates of thermal gradient than shown in Figure 5; we 
comment on this point further below. 

The flow law for crustal material is probably less well 
constrained in general than that for the mantle. To explore at 
least partially the effect of this uncertainty, we also consid- 
ered the flow law of anorthosite reported by Shelton [1981], 
Values of T m and dTIdz differ by about 5-15% from those 
obtained with the diabase flow law of Caristan [1982] for 
strain rates in the range JO -16 to 10 19 s 1 . The straight-line 
approximation to the ductile portion of the strength enve- 
lope, while a reasonable simplification for the olivine flow 
law [McNutt and Menard , 1982], with flow laws for crustal 
material leads to an underestimate of ductile strength in the 
lowermost mechanical lithosphere and thus a slight increase 
in estimated 7 m , for an assumed constant strain rate. 

A considerable uncertainty in the application of relations 
such as those depicted in Figures 5 and 6 arises from our 
poor knowledge of the thickness of the Martian crust. 
Long- wavelength topography and gravity, if intepreted in 
terms of a crust of uniform density and variable thickness, 
indicate a mean crustal thickness of at least 30 km [Bills and 
Ferrari , 1978], which corresponds to zero crustal thickness 
beneath the Hellas basin. Models of the Viking line of sight 
(LOS) Doppler radio tracking residuals over the Hellas basin 
and the 370-km-diameter crater Antoniadi are consistent 
with essentially local Airy compensation if the crust beneath 
these impact features is 120-130 km thick [Sjogren and 
Wimberly , 1981; Sjogren and Ritke, 1982]. From the models 
of Bills and Ferrari [1978] a 130-km-thick crust beneath the 
Hellas basin [Sjogren and Wimberly , 1981] corresponds to a 
globally averaged crustal thickness of about 150 km. LOS 
tracking data over Elysium Planitia and Olympus Mons can 
be fit with varying degrees of Airy isostatic compensation 
and crustal thicknesses of 30-150 km [Janie and Ropers , 
1983; Janie and Jannsen, 1986]. Available gravity data are 
thus permissive of a mean crustal thickness anywhere in the 
range 30-150 km, an interval consistent with, but not sub- 
stantially narrowed by, the predictions of differentiation 
models of an early Martian magma ocean [ Warren , 1988]. 
Whatever the mean crustal thickness, local variations of ±30 
km or more are also likely [Phillips et al ., 1973; Bills and 
Ferrari , 1978]. 

Given these uncertainties, our approach is to assume that 
the large values of elastic lithosphere thickness ( T e > 100 
km) determined from the local response to the Isidis mascon 
and Olympus Mons and from the global response to the 
Tharsis rise exceed the thickness of the Martian crust, while 
the values of 20-30 km obtained for T e beneath the Tharsis 
Montes and Alba Patera (Table 1) are less than or compara- 
ble to the thickness of the crust. This assumption is consis- 
tent with known constraints on the crustal thickness and, 
because of the lesser creep resistance of crustal material, 
results in minimizing the spread in inferred thermal gradi- 
ents. 

Finally, it should be noted that the estimation of mechan- 
ical lithosphere thickness and thermal gradient is made 
under the assumption that flexural stresses dominated the 
local lithospheric stress field at the time of load-induced 



Solomon and Head: Heat Flux and Internal Dynamics on Mars 


11,079 


[1985] and Laid et al. [1986] have inferred that the mantle 
plus crust of the SNC parent body has a U abundance of 16 
ppb a Th/U ratio of 3.5, and a K/U ratio of 2 x 10 For an 
inferred mass fraction of 22% for an Fe-FeS core [Drc,6^ 
and Wdnke, 1985; Laid et al., 1986] the U abunda ^ e 
expressed as a fraction of the whole planet is 12.5 pbb. The 
present heat production for such bulk abundances is equiv- 
alent, under steady state, to a mean heat flow of 14 mW m • 
Treiman et al. [1986] obtain element-element correlations 
Liter to those of Dreibus and Wdnke [1985] and Laulet al. 
[1986], with slightly smaller Th/U and KAJ ratios (3.0 and 
1 x 10 4 respectively), which for the same bulk U abundance 
would give a heat production about 20% lower. 

These values of heat flow may be compared with the 
thermal gradients in Table 1 and corresponding estimates ot 
thermal conductivity. For Tharsis Montes and Alba Patera 
T m is less than the crustal thickness, so we adopt a thermal 
conductivity appropriate for crustal material, 2. 

K' 1 [Clark, 1966]. Lithospheric thermal gradients of 10- 
Kkm' 1 (Table 1) correspond to heat flow values of 25-35 
mW m -2 . For the mechanical lithosphere supporting Olym- 
pus Mons, the Isidis mascon, and the Tharsis rise, most o 
the lithosphere probably consists of mantle matenal, so we 
adopt the thermal conductivity g ,v e« by Scja/z andS ^ 
mons [1972] for olivine (Fo 86 Fa 14 ) at 560-930 K. ™ 

K 1 . Mean lithospheric gradients less than 5-6 K 
(Table l) then correspond to heat flow values less than 20-24 
mW m' 2 . These last estimates of heat flow, and the associ- 
ated temperature distributions depicted in Figure 7, shou 
be reduced to be consistent with the likely contrast m 
thermal conductivity and gradients across the crust-mantle 
boundary, but in order to do so a crustal thickness mu 
assumed. For a crust 30-50 km thick and the above values of 

crustal and mantle thermal conductivity the upper bound on 

heat flow consistent with mean lithospheric gradients 1 
than 5-6 K km' 1 is 17-21 mW m . A global elastic 
lithosphere at least 100 km thick, indicated by flexural 
models for the support of the Tharsis rise [Willemann and 
Turcotte, 1982; Banerdt et al., 1982], implies lithospheric 
gradients of 7 K km' 1 or less. The corresponding upper 
bound on heat flow is 23-25 mW m' 2 for a mean crustal 
thickness of 30-50 km. 

Thermal gradients at or near the upper end of the range 
allowed for Olympus Mons, Isidis, and the globally averaged 
response to Tharsis thus give heat flow values similar to or 
slightly larger than that expected from U, Th, and K abun- 
dances in INC meteorites under the steady state assump- 
Uon Such upper bounds on thermal gradients are also 
JS a modest contribution to heat ow from 
secular cooling of the mantle [e.g., Sihu ere , 

Most published thermal history models for Mars, in contrast 
predict a heat flow larger than that permitted by these upp 
bounds, a result attributable to an overest.rnat.on of hea 
production in the interior. We suggest that the global heat 
flux on Mars during the Amazonian epoch (approximately 
the last 2 b.y. [Tanaka, 1986]) has been in the range l.-^ 
mW m' 2 This inference is in agreement both with the 
thermal gradients implied by the thick elastic llthos P her . e 
beneath Olympus Mons and globally supporting the Tharsis 
rise and with the bulk abundances of heat-producing el 
ments in Mars inferred from SNC meteorites. Regona 
flow in the centers of major volcanic provinces of Mars 


during the Amazonian epoch reached highs of 35 mW m 

more, as much as twice the typical global value. 

A final point of comparison with the thermal gradients ,m 
Table 1 and Figure 7 is provided by estimates of the 
thickness of the thermal lithosphere made from the heights 
SS constructs u„d a hydrostatic model for magma, tc 
overpressure [Vogt, 1974; Carr, 1976; Blastus and Cults 
1976] Such estimates are quite approximate, given 
tainties contributed by differential compress, b.l.ty between 
magma and surrounding rock, viscous head loss in the 
magma conduit, and additional overpressure contributed by 
magmatic volatiles, as well as by the probably oversimplified 
notion of a continuous magma column extending 
basal magma chamber to the volcano summit. We consider 
such models here merely as tests of consistency. For n 
stance the thermal gradient shown in Figure 7 for Olympus 
Mons, if continued downward, reaches a tem P e ‘‘ a ^ ire ° 
1670 K, the 2.3-GPa solidus temperature of a model Mar 
mantle [Bertka and Holloway, 1988], at a depth of 285 km. 
Inasmuch as the volcano stands about 24 km higher than e 
surrounding terrain [Wu et al, 1984] this relief As 
with the thermal gradient shown and with the hy 
head model if the average effective density -^^"Te 
maerna and the surrounding rock column is about 8%, a quite 
reasonable value. The lesser relief (14-18 km) of the Tharsis 
Montes volcanoes [Blast, is and Cults 1976] is _consis en 
with the higher thermal gradients indicated for thelitho 
sphere beneath these features (Table I and Figure 7), b» 
quantitative agreement with the magma hydi-oslattc model 
would require consideration of distinct densities and therm 
conductivities for the crust and 

crustal thickness [Blasius and Cutts, 1976], and ‘ * e P ° 
contribution of convective flow to heat transport in the lower 
lithosphere beneath the central Tharsis region. 


Evolution of the Tharsis Province 

As noted above, considerable attention has been devoted 
totheTep structure and evolution of the Tharsis volcanic 
and tectonic province. While the long-wavelength gravity 
and topography of the region are not consistent with cona- 
plete local isostatic compensation by a single mechams 
slich as crustal thickness variations [Phillips and ganders 
1975] complete local compensation is possible if a com 
Tuition of Airy (crustal thickness variations, ^ 
density variations) mechanisms act in concert, but only if the 
crust fs relatively thin (or is pervasively intruded by high- 
density plutonic material) beneath the Tharsis ^^an 
substantial density anomalies persist to at least 3<XM00 km 
depth [Sleep and Phillips, 1979, 1985; Finnerty et al., 1988]. 
Alternatively, the gravity and topography are consistent 
with the hypothesis that a portion of the high topograp y o 
Tharsis is supported by membrane stresses in the Martian 
elastic lithosphere [Willemann and Turcotte, 1982, B 
et al., 1982; Sleep and Phillips, 1985]. 

These compensation models have been used to ca cu a 
lithospheric stresses for comparison with the observe « 
tribution of tectonic features. The isostat.c models for Tta* 
sis predict stresses in approximate agreemen 
distribution and orientation of extensional fractures in th 
central Tbarsi, region and of compress, ve Wrinkle ndg es 
oriented approximately circumferential to the center o 
tectonic activity, while the models involving lithospheric 
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modest areas immediately surrounding localized lithospheric 
loads. The areas within major volcanic provinces for which 
the enhanced thermal gradients of Table 1 and Figure 7 are 
appropriate, unfortunately, are poorly constrained None- 
theless, an approximate measure of the lateral extent of 
enhanced lithospheric heating is the area of regionally ele- 
vated topography arising either from remaining excess hea 
or more likely, from a crust permanently thickened by the 
addition of volcanic and intrusive material. We assume that 
the areas of excess thermal gradient in the central Tharsis, 
central Elysium, and Alba Patera regions are 3 x 10 km , 

0.9 x 10 6 km 2 , and 0.8 x 10 6 km 2 , respectively. These areas 
correspond to the regions encompassed by the 7.5-km, 1- m, 
and 5-km elevation contours in the respective provinces 
[U.S. Geological Survey , 1989], The 7.5-km contour in the 
Tharsis region encloses the area of the Tharsis Montes, u 
not Olympus Mons, for which a lithosphere thickness and 
thermal gradient nearer to the global average values are 
indicated (Table 1). For excess thermal gradients of about 10 
K km -1 in central Tharsis and 5 K km m the other two 
volcanic provinces (Table 1) and a crustal thermal con uc- 
tivity of 2 5 W m' 1 K" 1 the heat delivered by these excess 
gradients is a, , raie of I x ,0" W, wUh 80% of the 
contribution to this figure coming from the Tharsis region. 
This figure is 3-5% of the estimated present global heat loss, 
2-4 x 10 12 W for the probable range of average heat flux 
values given above. By way of comparison the fraction of 
mantle heat flux delivered by plumes at present on Earth has 
also been estimated to be about 5% [Davies, 1988, Sleep, 
19901 The heat delivered by plumes on Mars must be 
supplied by cooling of the Martian core. In the thermal 
history models of Stevenson el al. [1983] the Martian core 
delivers heat to the mantle at a rate of 2 x 10 W during the 
Amazonian, sufficient to supply the estimated heat trans- 
ported by plumes beneath major volcanic provinces. 

We may compare these figures to the fractional heat flow 
delivered by volcanism and associated igneous intrusions. 
Estimates of the surface area of volcanic material at each 
major stratigraphic stage, including corrections for later 
burial, have been given by Greeley [1987] and Tanaka et at. 
[1988] Both find 2 x 10 8 km 2 of volcanic material, though 
the two analyses differ in detail, particularly in the relative 
strength of a “peak” in the flux curve at early Hesperian 
times (corresponding to the formation of the Martian ridged 
plains) about 3-3.5 b.y. ago [Tanaka, 1986]. Greeley [1987] 
has suggested that the volume of volcanic material may be 
estimated by multiplying the area by an average thickness of 
about 1 km. A volume V of 2 x 10 8 km 3 may be converted 
to equivalent heat Q by means of the relation Q - p(t p M 
A H)V where p is the density of volcanic material, C p is e 
specific heat, AT is the difference between eruption and 
ambient temperatures, and A H is the heat of fusion. We take 
p = 3 Mg m" 3 , C p = 1.2 kJ kg' 1 K-', and A H= 0.4 MJ 
kg" 1 , and we adopt AT = 1450 K [Bertka and Holloway, 
1988]. Averaged over 3.8 b.y., 2 x 10 km o ^vocamc 
material delivers only 10 10 W, or less than f ’ ° e 
estimated present global heat loss of 2-4 x 10 W. Accom- 
panying any volcanic eruption, of course is usually a 
significant intrusion of magma that cools and solidifies at 
depth The ratio of intruded to extruded volumes is as great 
as 10:1 in terrestrial eruptions [Crisp, 1984], The combined 
volumes of volcanic and intrusive material on Mars, for this 
ratio, would have delivered heat at an average rate of 3-6% 


of the total global heat loss. Note that the intrusive compo- 
nent of igneous activity may find expression in the inferred 
conductive gradient, depending on the depth of intrusion, 
but the volcanic component of activity will not. 

The volcanic flux has decreased sharply with time over 
Martian history [Greeley, 1987; Tanaka et al.,\9S ], so it is 
worth making a similar calculation both for the periods of 
high volcanic output as well as for the Amazonian epoch 
corresponding to most of the estimates of lithosphere thick- 
ness and thermal gradient given in Table 1 and Figure 7. Fo 
both the mid-Noachian and early Hesperian epochs, volca- 
nism resurfaced large areas in widespread regions over the 
planet. With the areas of volcanic material given by Tanaka 
et al [1988] for these periods, the 1-km thickness estimate o 
Greeley [1987], and the time intervals given by Tanaka 
[1986] 'the heat delivered by volcanism was at rates of 2 x 
10 n and 3 x 10 10 W, respectively. For a 10:1 ratio ° 
intrusive to erupted material the heat flux delivered by 
magmatism during the mid-Noachian may have been as 
much as one third of the global heat flux, with the precise 
ratio depending on mantle heat production, global thermal 
history, and the uncertain duration of the epoch. In contrast, 
the heat delivered by volcanism and plutonism estimated by 
this same procedure for the Amazonian contributed on 
average no more than 4 x 10 10 W, a figure equal to about 
40% of the excess conducted heat at major volcanic centers, 
estimated above, and only 1-2% of the global heat loss. 


Conclusions 

Estimates of the thickness of the effective elastic litho- 
sphere of Mars appropriate to a variety of locations an 
times have been converted to estimates of lithospher e 
thermal gradient and surface heat flow by means of strengt 
envelope considerations. Local thermal gradients and heat 
flow values were 10-14 K km and 25 35 mW m a 

time of formation of load-induced graben surrounding the 
Tharsis Montes and Alba Patera, while gradients and heat 
flow values of less than 5-6 K km - and 17-24 mW m 
characterized the lithosphere beneath the Isidis mascon and 
Olympus Mons at the time of emplacement of these loads. 
On the basis of the thickness of the global elastic lithosphere 
required to support the Tharsis rise, inferred on thermal 
grounds to characterize a late, rather than an early, stage 
the evolution of the Tharsis province [cf. Sleep and Philips, 
1985], as well as heat production estimates obtained from 
SNC meteorites, we suggest that the present global ea ux 
on Mars is in the range 15-25 mW m" 2 . Approximately 3-5% 
of this heat flux during the Amazonian epoch has been 
contributed by excess conducted heat in the central regions 
of major volcanic provinces, a figure at least a fac or 0 
greater than the heat transported solely by volcanism and 
shallow igneous intrusions. This excess heat flux is plausibly 
attributed to the action of mantle plumes on the base of th 
lithosphere beneath volcanic province centers. The frac- 
tional mantle heat transport contributed by plumes during 
the last 2 b.y. on Mars is therefore comparable to the present 
situation on Earth. 
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ABSTRACT 

The formation of a large volcano loads the underlying lithospheric plate and can lead to 
lithospheric flexure and faulting. In turn, lithospheric stresses affect the stress field beneath and 
within the volcanic edifice and can influence magma transport Modeling the interaction of these 
processes is crucial to an understanding of the history of eruption characteristics and tectonic 
deformation of large volcanoes. We develop preliminary models of time-dependent stress and 
deformation for the Tharsis volcanoes on Mars. A finite element code is used that simulates 
viscoelastic flow in the mantle and plate buoyancy forces. We show the manner in which time- 
dependent stresses induced by lithospheric plate flexure beneath the volcanic load may affect 
eruption histories, and we suggest that the stress field may be relatable to tectonic features on and 
surrounding Martian volcanoes. After an initial load, flexurally-induced stresses grow with time 
and the principal stress directions in the volcano rotate as flexure proceeds. Magma is expected to 
propagate perpendicular to the least compressive stress axis. As a result of flexure, this axis 
rotates from a horizontal orientation to a nearly vertical one; thus magma propagation paths will 
tend to rotate from vertical to horizontal orientations. We suggest that at the later stages of flexure, 
this effect would tend to favor eruption sites on to the flanks of the volcano rather than the summit. 
Such a scenario is consistent with the photogeologically determined evolution of the Tharsis 
Montes. During flexure there are three regions where stresses become sufficiently large to cause 
failure by faulting (according to the Mohr-Coulomb criterion): at the surface of the plate just 
outward of the edge of the volcano, on the volcano's flanks, and inside the plate beneath the center 
of the volcano. Normal faulting is the dominant mode of failure predicted for the first region, 
consistent with circumferential graben observed around the Tharsis Montes and with the scarp at 
the base of Olympus Mons, interpreted as a large-offset, listric normal fault Failure in the second 
region is predicted to consist of thrust faulting, oriented mostly circumferentially on the upper 
flanks with a narrow annulus of radially oriented thrust faulting about midway downslope. 
Concentric terraces, interpreted by some workers as thrust faults, on the upper flanks of Olympus 
Mons may correspond to the predicted circumferential thrust features. Normal faulting, mostly 
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radially oriented, is predicted for the mid-plate zone of failure beneath the volcano. Under the 
premise that failure in this zone may be strongly influenced by pre-existing weak zones or regional 
stresses, this feature may have a surface expression in the rifting and the bilateral symmetry of the 
Tharsis Montes. 
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INTRODUCTION 

The major Th arsis volcanoes of Mars (Olympus Mons and the three Th arsis Montes: Arsia 
Mons, Pavonis Mons and Ascraeus Mons) are among the largest known volcanic structures in the 
solar system. An understanding of the formation and evolution of these structures can provide 
important constraints on the processes that built and maintained the Tharsis rise. In addition, the 
Th arsis volcanoes may be analogues to large hot-spot volcanoes on Earth, such as Kilauea on 
Hawaii. For instance, Kilauea and Olympus Mons have very similar ratios of volcano height to 
diameter and of basal scarp height to volcano height ( Borgia etal., 1990). Thus studies of the 
evolution of large Martian volcanoes may yield insight into terrestrial volcanic processes as well. 

In this paper we utilize finite element models to evaluate the evolution in internal stress and 
deformation within and surrounding the Tharsis volcanoes, and we discuss how the time- 
dependent stress field may be relatable to the eruption characteristics of the volcanoes and to the 
formation of associated tectonic features. 

To date the investigation of the evolution of stresses in large volcanoes has taken two paths: 
models of edifice stresses alone (usually finite element models with rigid bottom boundary 
conditions), and investigations of flexural stresses in the lithospheric plate supporting the volcano. 
As examples of the first category of study, Chevallier and Verwoerd (1988) used an axisymmetric 
planar finite element code to investigate the effect of magma chamber and external pressures on 
stress and eruption histories of hot spot volcanoes, Dieterich (1988) modeled stress in volcano rift 
zones by means of a two-dimensional triangular grid of elements, and Ryan (1988) employed a 
horizontal planar finite element model of the flank of Kilauea to determine displacements and 
stresses due to dike emplacement In an example of the second class of study, Thurber and Gripp 
(1988) applied a flexure model to constrain the tectonics of volcano flank motions, ten Brink and 
Brocher (1987) proposed a link between flexural stresses in the lithosphere and eruption history; in 
their scenario the orientation of flexural stresses- beneath a given point on the volcanic chain 
changes with time as the volcanic load is emplaced and then eroded. Given that magma seeks to 
propagate along paths perpendicular to the least compressive stress direction, magma ascent can be 
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blocked when both principal horizontal stress deviators arc compressional. Thus eruption history 
depends on location along the volcanic chain or on time since volcano formation. 

Recently studies have been carried out on edifice stresses of large Martian volcanoes. Zuber 
and Mouginis-Mark (1990) treated the evolution of the surface of Olympus Mons, using a finite 
element model to calculate stresses in the summit caldera region caused by different magma 
chamber locations and geometries and to compare those stresses with patterns of faulting. Thomas 
et al. ( 1990) investigated the tectonics of the flanks of Martian volcanoes by means of an 
incompressible finite-element model. Volcano self-loading and magma chamber effects were 
included, but lithospheric flexure was not. Thomas and coworkers found that stresses on the 
flanks of a large volcano are sufficient to cause circumferentially oriented thrust faulting, and they 
suggested that such thrusting produced the concentric terraces observable high on the flanks of 
Olympus Mons. 

In this paper we study the stress field within a volcano and the lithosphere upon which it rests 
as a unified system. With such a model formulation we can account explicitly for the interaction 
between edifice stresses (e.g., due to volcano self-loading) and flexural stresses (a result of the 
load induced by the volcano on the lithosphere). We include viscoelastic deformation in the 
asthenosphere, so the problem is intrinsically time-dependent The calculated displacements yield 
the subsidence history of the volcano. The orientations of principal stresses and their change with 
time provide important constraints on possible magma emplacement paths and eruption histories. 
With the computed stress fields, a failure criterion can be used to predict locations and modes of 
faulting within and near the volcano. After a short discussion of important geological and 
geophysical characteristics of the Tharsis volcanoes, we briefly describe the finite element 
procedure and the modelling assumptions. The results of the numerical computations are next 
presented, and their potential implications for the evolution of eruption characteristics and for the 
formation of tectonic features are compared with known constraints on the evolution of the Tharsis 
volcanoes. 
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CHARACTERISTICS OF MAJOR THARSIS VOLCANOES 

The four major Tharsis volcanoes are the largest of the martian shield volcanoes (e.g., Greeley 
and Spudis, 1981). Each construct is composed of many overlapping flows and flow units 
erupted over a period of activity as much as 2-3 b.y. in duration ( Tanaka , 1986). The approximate 
heights and basal diameters are given in Table 1. 

The largest of the major Tharsis volcanoes is Olympus Mons (Figure 1). The principal 
tectonic features associated with Olympus Mons are (1) the summit caldera complex, consisting of 
a series of circular depressions with complex patterns of faulting, (2) the basal scarp, a nearly 
vertical cliff surrounding the volcano at a radius of about 300 km. Of less certain origin are (3) the 
aureole deposits, which occupy a vast region dominandy to the northwest of the main shield 
(downslope from the Tharsis rise), and (4) concentric terraces seen on the upper slopes of the 
volcano. The summit caldera is likely the result of a collapse following withdrawal of magma from 
a high-level magma chamber {Mouginis-Mark, 1981; Zuber and Mouginis-Mark, 1990). The basal 
scarp has been variously interpreted as a thrust fault ( Morris , 1981), a listric normal fault ( Francis 
and Wadgc, 1983), and a fault-propagation fold over a subsurface thrust fault ( Borgia et al., 

1990). The aureole deposits are generally held to be disrupted landslide material derived from the 
slopes of the volcano ( Harris , 1977; Lopes et al ., 1980; Francis and Wadge, 1983). As noted 
earlier Thomas et al. ( 1990) have interpreted the concentric terraces to be thrust faults. 

Tectonic features observed on and around the Tharsis Montes volcanoes differ somewhat from 
those of Olympus Mons. Noteworthy are the circumferential graben which occur on the lower 
volcano slopes and the surrounding plains (Figures 2-5). The Tharsis Montes volcanoes exhibit a 
qualitative bilateral symmetry about a NE-SW-trending axis coinciding approximately with the line 
connecting their centers. From photogeological study of Viking Orbiter images, Crumpler and 
Aubele (1978) proposed the following evolutionary sequence for the Tharsis Montes: (1) 
construction of the main shield, (2) outbreak of parasitic eruption centers on the volcano flanks 
along the NE-SW-trending axis, (3) subsidence of the summit and formation of concentric 
fractures and graben, and (4) formation of a bisecting rift along the NE-SW-trending axis, with rift 
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eruptions leading to flooding of the summit depression and inundation of the rifted flanks. This 
evolution is most developed on Arsia Mons; Pavonis Mons has reached stage (3), and Ascraeus 
Mons stage (2). 

Comer et al. (1985) inferred the thickness of the elastic lithosphere beneath the Tharsis 
Montes volcanoes from the radial distances of their circumferential graben (see Figures 3-5). 
Preferred values of elastic lithosphere thickness are around 20 km. Concentric graben are not 
found around Olympus Mons, which led Comer and others to conclude that the elastic lithosphere 
thickness under Olympus Mons must be much greater (> 150 km) than beneath the Tharsis 
Montes. Of course, the actual lithosphere does not behave perfectly elastically; rather its strength is 
limited by frictional failure at shallow depth and by ductile flow at greater depth. From strength 
envelope considerations for crustal and mantle material ( McNutt , 1984), values for elastic plate 
thickness were converted into estimates of mechanical plate thickness by Solomon and Head 
(1990). 


METHOD 

We use the finite element program TECTON ( Melosh and Rqfesky, 1980, 1983) to model 
stresses and displacements in a large volcano and in the crust and mantle beneath and around the 
volcano. TECTON's capability for adopting a viscoelastic rheology allows us to model time- 
dependent plate flexure. The program first calculates the elastic (i.e., "instantaneous") response to 
the load. The stresses and displacements arising from load-induced viscoelastic flow in the mantle 
are then calculated for a specified number of time steps. The Maxwell time (x m) of the mantle, 
defined as the ratio of viscosity t| to shear modulus |i, is used as a convenient reference time scale. 
For all computations here, timesteps were implemented in three groups. The first group had 
timesteps of duration about equal to x m and were run up to 10 x the second group had steps of 
10 x m run up to 100 x m , and the third group used steps of 100 x m ending at 1000 x*m . Results 
shown here were obtained at approximately 10, 100, and 1000 Maxwell times. 
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We assume that the problem is axisymmetric, with cylindrical coordinates r, 0, and z. Out of 
plane shear stresses G^ and Gq z are then zero. We solve for r and z displacements and stresses 
G rr , Gqq, g zz , and G ri . In axial symmetry, two principal stresses are confined to the rz plane; thus 
CTq 0 must also be a principal stress. An example of the finite element grid used for this study is 
shown in Figure 6 . The displacement boundary conditions are that nodes on the side walls (r = 0 
and r = r max ) are fixed in r but free to move in z, and that nodes on the bottom boundary (z = -940 
km in the example shown) are fixed in z but free to move in r. The lower comers of the box are 
fixed in both directions. The volcano in the example has a radius of 200 km and is 20 km in 
height, the approximate present dimensions of Ascraeus Mons. In Figure 6 the volcano is the 
triangular region in the upper left-hand comer. The volcano rests on top of an elastic lithospheric 
plate of thickness T e . All elements in the volcanic edifice and the plate have a high viscosity 
appropriate to the lithosphere. These elements behave essentially elastically over the timescales 
considered here. All elements below the elastic lithosphere have a lower viscosity value, 
appropriate for an asthenosphere. These elements experience viscous relaxation over time. 
Material property values adopted in our calculations are listed in Table 2. Parameters such as 
density and Young's modulus differ for the crust and mantle. Elements above depth ^ arc 
assigned crustal values for these parameters; below this depth mantle values are assigned. 
Somewhat arbitrarily, for models discussed here, we have chosen T c = tj.. We do not mean to 
imply that thicknesses of the Martian lithosphere and crust coincide; this choice serves only to 
simplify the model. The results of three computations are discussed here, with T e values of 20, 

40, and 80 km. These values chosen for T e fall within the range discussed in Comer et al. (1985) 
for the Tharsis Montes and Olympus Mons. For T e = 20 and 40 km, the outer radial boundary of 
the grid is taken as r max =1140 km. For T e = 80 km, such a value was found to be too small to 
avoid end effects, so we took r max = 1600 km. 

In order to model plate flexure correctly, we must include the effect of buoyant support of the 
crust by the denser material beneath it In TECTON, this is accomplished through the use of 
Winkler restoring forces at the nodes marking the boundary between the crust and the mantle (at a 
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depth t(.). These Winkler forces are proportional to displacement. The buoyant force (per unit 
area) exerted by the mantle on the plate is: 

> 

fb = 'Pm 8 5z 

where p m is the mantle density, g is the gravitational acceleration (3.7 m/s 2 ), and 5 z is the vertical 
displacement. We specify Winkler forces of this magnitude at each node along the crust-mantle 
interface (applied over the appropriate area for that particular node). 

The calculations performed here should be regarded as simple numerical experiments designed 
to explore the stress field resulting from a volcano-like load on an elastic plate underlain by a 
viscoelastic substrate. When applying the results of these calculations to actual geologic and 
tectonic features seen on Martian volcanoes, one must keep in mind the limitations of such models. 
The first and perhaps most important limitation is the instantaneous placement, at time t = 0, of a 
fully developed volcanic edifice. This simplification implies a volcano growth time t g much less 
than the mantle Maxwell time x^. For Mars, if T| ~ 10 21 Pa-s and p. ~ 10 1 1 Pa, x^ — 10^ s or 
300 yr. The condition t g « is thus clearly unrealistic. The effect of episodic increments to 
load size (due to episodic eruptions) over many Maxwell times are not amenable to investigation 
with this approach. In addition, even for the assumption of an instantaneous load, it is difficult to 
determine appropriate values for the initial geometry of a given volcano, since existing data apply 
only to volcanoes that have already undergone flexure and deformation that may have significantly 
modified initial conditions. These calculations are performed under the assumption of axial 
symmetry; possible effects of nonaxisymmetric loading (such as regional stress) and complex 
three-dimensional geometry are not incorporated. The effects of magma pressure, transport, and 
evacuation are not addressed in these preliminary models; these processes likely have important 
influences for caldera and flank tectonic evolution ( Thomas et al., 1990; Zuber and Mouginis- 
Mark , 1990). Further, lateral variations in material properties (due to horizontal temperature 
gradients, for example) are not included. Our choice that T e = t,., while made for convenience, is 
not necessarily valid for the Tharsis region. Finally, it must be remembered that once faulting 
occurs, the stress fields calculated are no longer strictly valid, since faulting would relieve stresses 
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locally and could thus alter the predictions for failure at subsequent times. As discussed below, 
some of these limitations can be relaxed in future modeling efforts. 

NUMERICAL RESULTS: STRESS, DEFORMATION, AND FAULTING 

The evolution of the stress field within and beneath the volcanic edifice in a model is dipicted 
by means of symbols for the principal stresses within each element (Figure 7). The hourglass 
shapes are oriented along directions of deviatoric compression, while the bars are oriented along 
directions of deviatoric tension; symbol dimension is proportional to stress magnitude. For a 
general state of stress in a plane layer (without the volcanic load), the orientation of the stress 
symbol will be as in Figure 7: maximum compression vertical and maximum extension horizontal. 
The magnitude of the deviatoric stress will increase with depth. For a medium constrained 
laterally, this result can be derived from the equations for uniaxial strain (compaction) in the z- 
direction (e.g., Turcotte and Schubert , 1982 , p.108): 

°rr * c 98 = (jt^) a zz 

For v = 0.25, the factor v/(l - v) is 1/3. Thus, the vertical stress will be three times as great as the 
horizontal stress, and with vertical stress equal to the overburden stress (G zz - p c gz\ the stress 
difference will increase with depth. 

One such model for the evolution of volcano-related stresses is shown in Figure 8. At t = 0 
(Figure 8a) deviatoric stresses display the orientations and magnitudes expected from simple self- 
loading of horizontal layers everywhere but along the top layer of volcano elements, where the 
stresses are rotated such that the compressive axis is almost horizontal. As the effects of flexure 
manifest themselves, this rotation propagates deeper into the volcano, eventually reaching the 
upper part of the underlying plate (Figure 8b-d). 

m 

The evolution of surface elevation, relative to its value at the right boundary r = r max , is 
shown as a function of time in Figure 9, for computations with three different values of T e . For all 
cases, most of the subsidence takes place between 10 and 100 Maxwell times. The maximum 
amount of subsidence occurs over the center of the volcanic load. This subsidence results in a 



11 


reduction of volcano slope with time. Table 3 lists the average volcano slope for each of our 
computations for the elastic solution and as a function of time. The slope for the assumed initial 
volcano has a value of 6.34°. The amounts of subsidence and slope change increase with 
decreasing plate thickness. 

The magnitudes of the surface stresses G^, and Gqq are shown as functions of time in Figure 
10. The stresses depicted are actually calculated at the center of each linear quadrilateral element at 
the surface of the model. The jaggedness in the curves occurs because the depth of the top 
elements changes along the slope of the volcano. 

Given stress values, we may use the Mohr-Coulomb failure criterion to estimate regions 
where faulting would occur. The Mohr-Coulomb failure equation relates shear stress X at failure to 
normal stress G n 

^failure = C + Gn tan 0 

where c is the cohesive strength of the rock and <j> is the angle of internal friction of the rock. We 
adopt values for c (3.8xl0 7 Pa) and <(> (49°) appropriate for basalt (Handin, 1966). Once we have 
found where failure is expected, the orientations of principal stresses are used to determine the 
style and orientation of faulting, according to the criteria set forth by Anderson (1951). Given the 
directions of the horizontal and vertical normal stress components, the type of faulting (normal, 
thrust, or strike-slip) and the orientation (radial or circumferential) can be determined. Table 4 
gives a listing of the possible types and orientations of faulting resulting from given stresses. Here 
Gj refers to the greatest compressive stress, 03 to the greatest extensional stress. Elements where 
failure predicted to have occurred are assigned a symbol corresponding to the stress state at that 
node. This classification scheme holds if the principal stresses approximately correspond to the 
horizontal (g^ Gqq ) and vertical stresses (g zz ) listed here. If the stresses are not aligned in this 
way, the asterisk symbol is assigned. Strike-slip faults do not quite fit the usual definitions of 
radial and circumferential, because such faults strike obliquely to the principal stress directions. 

For most rocks, strike-slip faults make smaller angles with the Gj than the G3 direction. Thus, 
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when a! = (Iqq we will consider the resulting faults circumferentially oriented, and when Oj = o rr 
we consider them radial. Care must be taken when applying these results to regions at depth, or 
where principal stresses are rotated obliquely to the surface. For the purposes of applying the 
categorization of Table 4, principal stresses making angles of less than 20° with the horizontal are 
considered horizontal, and similarly for near vertical stresses. 

The predicted locations and types of faulting as functions of time are shown in Figures 11-13 
for the three values of elastic lithosphere thickness. Figure 1 1 depicts the case with T e = 40 km. 
After 10 Maxwell times (Figure 1 la), a region of failure has initiated near the outer edge of the 
volcano surface. This region exhibits mostly circumferential normal faulting, with a small section 
of circumferential strike-slip faulting toward the volcano. After 100 Maxwell times (Figure 1 lb), 
this region expands dramatically. Two new zones of failure also appear a zone of thrust faulting 
along the upper flank, and a region of normal faulting of an average depth of 30 km beneath the 
center of the volcanic load. The predicted fault orientation on the flanks is predominantly 
circumferential, but a radial orientation is predicted for a small annulus approximately midway 
downslope. The region of subsurface faulting lies between 20 and 40 km depth and extends 
laterally beyond the outer volcano radius. It consists of a core region of circumferential normal 
faulting, surrounded by a region of radial normal faulting. After 1000 Maxwell times (Figure 
1 lc), little change is evident other than a slight outward extension of the surficial region of normal 
faulting and a greater volume of mixed fault geometry in the subsurface faulting region. 

Changing the plate thickness can affect the evolution of failure regions. For T e = 20 km, 
failure is predicted in only one element after 10 x M (Figure 12a). After 100 x M (Figure 12b), 
regions of faulting similar in location and style to those in Figure 1 lb appear, with the 
modifications that the region of strike-slip faulting is shifted inward and up the volcano flank, the 
region of surficial normal faulting is smaller in extent, and the region of subsurface normal faulting 
is shallower (10-20 km deep). Also, the subsurface normal faulting region takes longer to extend 
beyond the volcano radius (Figure 12b). 
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Results for the computation with T e = 80 km also reveal interesting changes. At t = 10 
(Figure 13a), the only faulting occurs at the surface, at the edge of the volcano. At t = 100 
(Figure 13b), flank thrust faulting is predicted, similar to that in Figure lib. The subsurface 
faulting region is found at a deeper level (60-80 km depth ). The outer faulting region shows the 
most dramatic change from the calculations with thinner lithospheres. This region expands 
outwards to more than two volcano radii from the center, as well as down to a maximum depth of 
20 km. Circumferential strike-slip faulting is much more prominent than in the previous cases, and 
is farther from the volcano. Finally, at t = 1000 (Figure 13c), we observe that the outer region 
of faulting has migrated even farther, out to a maximum extent of three volcano radii from the 
center, and that the area undergoing strike-slip faulting has grown slightly. 

DISCUSSION: IMPLICATIONS FOR VOLCANO EVOLUTION 

The volcanic edifice stress fields calculated here have important implications for magma 
transport and eruption. Magma propagates through fractures that form perpendicular to the 
direction of least compressive stress. The elastic solution (Figure 8a) clearly shows that magma 
can propagate vertically, with a clear path to the summit caldera. As a result of flexure, however, 
the principal stress directions in the edifice rotate such that the most compressive stress is 
horizontal. In this region, the magma propagation direction should tend to the horizontal, leading 
outward from the summit to the flanks. Thus an implication of the models treated here is that there 
should be an evolution in favored eruption location from summit to flank. Such an evolution is 
consistent with the first two stages in the sequence of major events for the Tharsis Montes 
determined by Crumpler andAubele (1978). 

The predicted types and orientations of faulting in our models (Figures 1 1-13) can be 
compared with observed features on and around the Tharsis volcanoes. Circumferential graben are 
observed on the lower flanks and the surrounding plains for all three of the Tharsis Montes 
(Figures 2-5). An interpretation of the basal scarp of Olympus Mons as a large-offset listric 
normal fault ( Francis and Wadge, 1983) is also consistent with model results. Circumferential 
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terraces high on the flanks of Olympus have been interpreted as concentric thrust faults by Thomas 
etal. (1990), who applied elastic finite element models to investigate edifice stresses and 
deformation. Interestingly, our models indicate that these faults would form only as a result of 
flexural stress, whereas the models of Thomas and coworkers do not include flexural stresses, but 
suggest that purely elastic stresses are sufficient to exceed the compressive failure strength at near 
surface conditions. The thrust faulting region in our model is somewhat difficult to reconcile with 
the observation of graben rather high on the slopes of Arsia and Pavonis Mons (Figures 4, 5). 
Clearly, an extensional environment existed in the upper volcano flanks at some time in the 
evolution of these constructs. These graben may be related, however, to faulting beneath the 
volcano, as discussed below. 

Linear symmetries observed in all the Tharsis Montes volcanoes may be indirect evidence for 
failure in the sub-volcano region. The second stage in the Crumpler and Aubele (1978) sequence 
of volcano evolution involves the development of a linear rift bisecting the volcanoes, with 
eruptions emanating from this rift Radially oriented normal faulting under a volcano might be 
thought at first to tend to divide the volcano into radial sections, much like pie slices. The above 
evidence, however, suggests that this mode of failure may be concentrated into one large linear 
feature that bisects the volcano. This feature may owe its orientation to regional stress fields or to 
pre-existing zones of weakness. Our calculations predict radially oriented normal faulting in a 
broad zone beneath a volcano. For the case with T e = 20 km (Figure 12), this zone extends 

beyond the volcano radius and reaches quite close to the surface. Such a zone, if concentrated into 
a linear feature and if continued to the surface by sustained faulting, may give rise to the linear rifts 
and bilateral symmetry observed on the Tharsis Montes. 

The evolution of Pacific hotspot volcanoes may provide some insight into the formation and 
influence of a bisecting rift zone. For instance, Holcomb (1985) has suggested that the island of 
Molokai in the Hawaiian chain resembles a volcano cut almost exactly in half, with one part 
transported downslope in a massive submarine slide, and the other part remaining to form the 
present island. Bathymetry of the ocean floor north of Molokai {Moore, 1964) supports this 


15 


conclusion. We conjecture that a similar but less extreme mechanism could explain the presence of 
graben on the upper slopes of Arsia Mons and Pavonis Mons, in that as the volcano halves slide 
apart from each other, an extensional environment is created near the summit To fully investigate 
this possibilty, a fully three dimensional model, or an axi symmetric model with non-axisymmetric 
loading, is required. 

We should also note that there are features predicted in our calculations that have no 
observable counterparts on the Tharsis volcanoes. In particular, evidence for strike-slip faulting 
around any of these structures is lacking. For values of lithosphere thickness appropriate to the 
Tharsis Montes, however, the zone of predicted strike-slip faulting lies at the edge of the volcanic 
construct, which has the lowest surface elevation (Figure 9). This topographic low could 
accumulate lava flows and erosional deposits which would cover the strike-slip zone, hiding it 
from photographic surveys. Alternatively, since the preferred initial failure mode for all except a 
very narrow annulus is concentric normal faulting, perhaps the initially formed faults continue to 
accommodate release of stress at subsequent stages in the development of flexure. Evidence for 
radially oriented thrust faulting has also not been observed on the lower flanks of any of the 
Tharsis volcanoes. Obscuration of such faults by flank lava flows is a possible explanation. If the 
rifts bisecting the Tharsis Montes are actually zones of radial normal faulting from which both 
halves of the volcano have slid apart, as suggested above, then the propagation of the normal 
faulting zone to the surface could change the stress fields in the lower flanks, relieving the 
compressive stresses expected from our calculations. 

The relative timing of events in our models is generally consistent with the chronology set 

forth by Crumpler and Aubele (1978) for the Tharsis Montes. After construction of the main 

► 

shield, rift formation and flank eruptions along this rift occur in the second stage of evolution. 
Figure 8b shows that, after 10 Maxwell times, the principal model result is the beginning of 
rotation of the stress axes in the edifice (which is expected to favor flank eruptions). After a few 
tens of Maxwell times the failure zone beneath the volcano (which could favor formation of the rift) 
has begun to grow. The third stage involves summit subsidence and formation of concentric 
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fractures and graben. After 100 Maxwell times, our computations show increased subsidence 
(Figure 9) as well as the growth of the region of normal faulting outwards of the volcano (Figure 
lib). The fourth stage involves further eruption along the rift, flooding the summit and flanks. 
Late stage flank eruptions are consistent with the stress fields in our calculations. Summit 
eruptions, however, would not be favored by the computed stress fields. It is likely that the stress 
field is influenced by the effects of previous stages of the evolution. Flank eruptions result in an 
added flank load, not modelled here, to which the lithospere will respond with additional flexure. 
Subsidence and additional flexure near the flanks of the volcano may at least partially relieve the 
compressional stresses near the summit, thus enhancing the possibility of additional summit 
eruptions. By this scenario, the preferred site for eruptions could alternate between summit and 
flank, on time scales on the order of 100-1000 Also, the effect of sections of the volcano 
sliding away from the bisecting rift may result in an extensional environment near the summit, once 
again favoring summit eruptions. 

Realizing the limitations of these calculations in accurately modeling the behavior of large 
volcanoes, we hope to be able to relax some of these limitations in future models. The restriction 
of small volcano growth time t g can be removed with incremental loading, a capability which 
requires code development Such a capability can also permit investigation of the effects of flexure 
on modulating eruption timing and location. A study of the effects of faulting on volcanic 
displacements and stresses can be accomplished through existing provisions in TECTON for 
modeling fault discontinuities, by means of the slippery node method ( Melosh and Williams, 

1989). 

CONCLUSIONS 

In this paper we have used finite element calculations that simulate the effects of lithospheric 
flexure to investigate the deformation and stress histories resulting from a volcano load on the 
lithosphere. The results of these calculations are compared with observed tectonic features and 
inferred eruption characteristics of large Martian volcanoes. Some of the zones of failure and 
faulting predicted in our models have direct analogues in observed features; others may be 
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indirectly related to actual structures. Three regions of failure with characteristic stress regimes are 
identified in our calculations. Stresses caused by flexure and subsidence may lead to thrust 
faulting on the flanks of large volcanoes. Flexural stresses in the lithospere immediately 
surrounding the volcano can result in circumferential graben. Flexurally induced failure in a wide 
zone beneath the center of a large volcano may play a major role in the development and 
modification of the edifice. Given that magma propagates perpendicular to the least compressive 
stress, a rotation in the principal stress orientations during lithospheric flexure in the models tends 
to favor eruption sites shifting from the volcano summit to the flanks as flexure proceeds. We 
conclude that time-dependent lithospheric flexure is important in determining the location and style 
of faulting within and surrounding large volcanoes as well as regulating the timing and location of 
volcanic eruptions. 
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FIGURE CAPTIONS 

Figure 1. Viking Orbiter view of Olympus Mons. The multiple pit caldera and basal scarp (at 
bottom) are the main tectonic features visible. No graben are seen on the flanks of the volcano 
or on the volcanic plains immediately surrounding the structure. Frame V0641 A52; width of 
image is 440 km. 

Figure 2. Viking orbiter view of western Ascraeus Mons. The complex caldera is seen at the 
summit, and graben with a predominantly circumferential orientation are seen along the 
western and southwestern margins. Graben are located on Ascraeus Mons itself and in 
surrounding flows (unit as, Figure 3). Frame V0643A78; width of image is 420 km. 

Figure 3. Geologic map of Ascraeus Mons and surroundings, simplified from Scott et al. (1981b) 
by Comer et al. (1985). Volcanic units shown include relatively young Ascraeus Mons flows 
(as), intermediate- age Tharsis Montes flows (tm), and volcanic material undivided by age (vu). 
Also shown as a distinct unit is slide material (s), interpreted by Scott and coworkers as 
landslides and debris flows. Dashed lines show approximate elevation contours, in kilometers, 
relative to a fourth-degree, fourth-order equipotential (Wu, 1978). The summit caldera 
complex is indicated by inward hatched lines. Extensional faults and graben are shown as 
heavy lines. 

Figure 4. Geologic map of Pavonis Mons and surroundings, simplified from Scott etal. (1981a, 
b, c) and Scott and Tanaka (1981) by Comer et al. (1985). Units shown, in addition to those 
described in Figure 3, are relatively young volcanic flows from Pavonis Mons (pm) and Arsia 

I 

Mons (am). Other information follows the format of Figure 3. Circumferential graben reach 
quite far up the slopes, and an approximate bilateral symmetry can be seen about an axis 
trending approximately NNE-SSW. 

Figure 5. Geologic map of Arsia Mons and surroundings, simplified form Scott et al. (1981c) by 
Comer et al. (1985). See Figures 3 and 4 for further explaination of symbols. Circumferential 
graben extend almost up to the summit caldera complex. 



22 


Figure 6. The finite element grid used for the calculations shown in later figures. The volcano is 
the small triangular area in upper left-hand comer. 

Figure 7. Key for stress symbols in Figure 8. The hourglass shapes are oriented along directions 
of deviatoric compression, the bars along directions of deviatoric extension. Symbol size is 
proportional to stress magnitude. Magma is expected to propagate perpendicular to the least 
compressive stress. 

Figure 8. Plots of deviatoric stresses in elements in the vicinity of the volcano, for the case with a 
40 km elastic lithospere thickness. The hourglass shapes denote deviatoric compression, the 
bars deviatoric extension, (a) The elastic solution, (b-d) After 10, 100, and 1000 Maxwell 
times, respectively. 

Figure 9. Plots of the elevation of the surface nodes of the model, relative to the surface node at 
the far end of the grid (r = (a) T e = 20 km. (b) T e = 40 km. (c) T e = 80 km. 

Figure 10. Plots of stresses calculated in the surface elements of the model, for the case with T e = 
40 km. Stress is taken to be positive in extension, (a) The elastic solution, (b-d) After 10, 
100, and 1000 Maxwell times, respectively. Note that the stress scale increases with each time 
increment. 

Figure 1 1. A plot of the elements at which failure is predicted, giving the style and orientation of 
expected faulting. The lithosphere thickness is T e = 40 km. (a-c) After 10, 100, and 1000 
Maxwell times, respectively. 

Figure 12. Plot of expected location and geometry of faulting versus time for T e = 20 km. (a-c) 
After 10, 100, and 1000 Maxwell times, respectively. 

Figure 13. Plot of expected location and geometry of faulting versus time for T e = 80 km. (a-c) 
After 10, 100, and 1000 Maxwell times, respectively. 



Table 1. Dimensions of Tharsis Volcanoes 


Volcano 

Diameter 

(km) 

Relief 

(km) 

Olympus Mons 

600* 

24 a 

Ascraeus Mons 

400 b 

18 b 

Pavonis Mons 

320 b 

l 4 b 

Arsia Mons 

420 b 

19 c 


a Wu etal. (1984). 
b Blasius and Cutts (1976). 
c Blasins and Cutts (1981). 


Table 2. Adopted Parameter Values 


Parameter 

Crust 

Mantle 

E, Pa 

lxlO 11 

3xlO n 

p, kg/m 3 

3000 

3500 

V 

0.25 

0.25 


Lithosphere 

Asthenosphere 

T|, Pa s 

lxlO 25 

1x1021 


Table 3. Average Slope Angles (in degrees) 


T e 

Elastic 

t=10 t m 

t=100 

t=1000 

20 km 

6.30 

6.10 I 

5.00 

4.32 j 

40 km i 

6.29 

~KU 1 

5.85 

5.80 

80 km 

6.28 

6.16 

5.85 

5.80 
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